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The Cluster Variation Method (CVM) is appHed to the Ishibashi model for ammonium di- 
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staggered and the uniform susceptibility without hysteresis are calculated at equilibrium. On 
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§1. Introduction 

Recently, the pyroclore oxide crystals A2B2O7 (A=Ho, Dy, B=Sn, Ti) called the spin ice 
have drawn many attentions of researchers. B The spin ice system has a typical geometrical frus- 
tration structure. The similar frustration structure is found in ammonium dihydrogen phosphate 
NH4H2PO4 (ADP). ADP is one of the hydrogen bonded crystals similar to well-known potassium 
dihydrogen phosphate KH2PO4 (KDP). However, as a result of the replacement of potassium ion 
by ammonium ion NH^, ADP undergoes the anti- ferroelectric phase transition while KDP is 
the typical ferroelectrics. Nagamiyai^ first suggested a possibility of anti-ferroelectricity for ADP in 
the framework of the Slater's model for KDP& in which the replacement of the first excited energy 
eq due to hydrogen configuration around PO4 by a negative value — eo might well explain the exper- 
imental behaviors in ADP. However, Ishibashi et o/.& pointed out that only taking negative value 
— eo is not enough to realize the antiferroelectric phase transition because there are several proton 
configurations with the same energy as that proposed by Nagamiya due to geometrical frustration. 
In order to single out the observed crystal structure in ADP experiments, Ishibashi introduced the 
dipole-dipole interaction in Nagamiya's proposed model. Ishibashii^ further analyzed the extended 
model (hereafter Ishibashi model) in which it is possible for three or four protons to come closer to 
a PO4 tetrahedron contrary to the ice rule. Here the ice rule demands: (1) each bond connecting 
the oxygen atoms in neighboring PO4 tetrahedra has always one and only one proton, and (2) 
each PO4 tetrahedron can have only two protons near-by. Namely the second rule of the ice rule 
is loosened in the Ishibashi model. The model with this type of extension for the Slater's KDP 
theory is often called the Slater- Takagi modelJl' 

In the present paper the Ishibashi model for ADP is reconsidered for the preliminary study 
of pyroclore crystals with geometrical frustration structure. We calculate the staggered and the 
uniform susceptibility above and below the transition temperature in the Ishibashi model for the 
ADP-type crystal in the cactus approximation of the cluster variation methodic which is equivalent 
to the Slater's approximation for the KDP model. Further, the hysteresis phenomena of uniform 
susceptibility versus temperature observed in experiments are successfully shown by making use of 
the natural iteration method (NIM) developed by KikuchiJl' As the case of hysteresis phenomena 
depending upon the external electric field, the polarization P curves against the uniform field is 
also calculated. 

In §2 we derive the variational free energy in the cactus approximation of the CVM.!-* In §3, from 
thermal equilibrium condition we obtain a self-consistent equation for polarization in a staggered 
electric field in order to find the properties of phase transition in the present system. After deter- 
mining the tricritical point which stands for the boundary between the first and the second order 
phase transition, we calculate the staggered susceptibility and the antiferroelectric sublattice spon- 
taneous polarization. In §4, after we calculate the uniform susceptibility of the present system in 
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thermal equilibrium, we study hysteresis phenomena of uniform susceptibility versus temperature 
in order to compare with the experimental results. §5 is devoted to a summary. 

§2. Formulation 

Let us consider the system consisting of 2N protons around PO4 tetrahedra in the ADP-type 
crystal as shown in Fig. 1. In order to formulate the present system in which the anti- ferroelectric 
phase transition along the a-axis takes place, we take a four sublattice model by Ishibashil^ as 
shown in Fig. 2. Here p, q, r and s denote non-equivalent hydrogen bonds on which a proton is 
located. When a proton on the hydrogen bond p is located close to the z-th sublattice, we denote 
the occurrence probability of such a proton configuration as pi. Next, let us assign the energy, 
the occurrence probability and the dipole moment along the a-axis to the protons configuration 
around PO4 tetrahedron as shown in Fig. 3. We apply the cactus approximation of the CVM to 
the present system so as to find the variational free energy. The cactus approximation of the CVM 
is equivalent to the Slater's theory for KDP which takes account of the site of a proton in the 
double well potential along each 0-0 bond (hydrogen bond) between two nearest neighbor PO4 
tetrahedra and the correlation of four protons around each PO4 tetrahedron. In the following we 
call the occurrence probability of a proton configuration on the hydrogen bond as a bond (state) 
variable and that of four protons configuration around PO4 tetrahedron as a four protons (state) 
variable. The number of configurations of L ensemble in the cactus approximation of the CVM is 
given byH* 

27V 

^=n<L n gs, (2.1) 

i=i {ijki) 

with 

TjAijkl) 

^(ijkl) _ tetra /r, r,\ 

^^bond ^^bond ^^bond ^^bond 

(i) 

where, for example, W^^J^^ is the number of proton configurations on the i-th bond if the i-th bond 
is p bond between 1 and 2 sublattice: 

^''-'^ = {Lp,y.{Lp2)V (^-^^ 

and cf^l^J"]^ is the correlation number of protons for a PO4 tetrahedron surrounded by protons on 

(ijkl) 

the i, j, k, 1 hydrogen bonds. And W^^l^gl is the number of four protons configuration around 
PO4 surrounded by protons on the i, j, k, 1 bonds. For example, referring to Fig. 3, for a PO4 

(ijkl) 

tetrahedron belonging to the 1-st sublattice, W^J^^J is given by 

[(L4^))!]^[(Lc«)!]^(La«)!(La«)![(La«)!]^[^^^^^^ ^ ' ^ 
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where the superscript (1) denotes the number of sublattice. By utiUzing StirUng's formula and 
assuming the homogeneity in each sublattice, the entropy of the present four sublattice model 
becomes 

S 1 

— = -logW 
Kb L 



N 

T 



^ \upi + Qi In Qi + Ti In n + Si In Sj 

i=l 

- (a« In + a« In a^^ + Ad^ In 
+4d« lnd«+2Gj') In 

+24'hn4')+24'hn4'))]. (2.5) 
Further, the electric polarization of i-th sublattice per PO4 along the a-axis is defined by 
/.„P« = - a«) + 2(d« - = 1 ~ 4), 

(2.6) 

where Ha is the dipole moment along the a-axis. The proton configuration energy U per system is 
given by 



N 
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1=1 

+ ^A/i2(p(l) + p(2))(p(3) + p(4)) (2.7) 

-^K^(p(i) + p(2)) + ^„£;'(p(3) + pW)], 

where the first line represents the protons configuration energy around PO4 tetrahedra, the second 
line denotes effectively the long range dipole-dipole interaction energy to induce an anti-ferroelectric 
structure along the a-axis and the last line is the energy due to external electric field. In the following 
we call the case of E' = —E a staggered electric field and the case of E' = E a homogeneous electric 
field. Here it should be noted that bond variables are always expressed in terms of four protons 
variables: 

2p, = 2q, = 4^) + 4') + a« + a« + ^d^ + d^l^ 

+4'^ + 4') + 4^) + + 3d? + 
2p, = 2q, = 4^) + 4^) + a(^) + a^ + d? + Sd^ 

+4'^+4'^+a?+aL^) + d?+3dL^), 
2n = 2.3 = 4'^ + 4'^ + + a? + 3d? + d« 

+4'^ + 4'^ + a?+aL^) + d? + 3dL^), 
2r3 = 2s, = 4^) + 4^) + a? + a? + d? + 3dL^) 
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2P3 = 2,4 = + cf + ai^) + aL^) + d^^) + SdL^^ 

27^4 = 2,3 = c^'^ + ^ + 0^ + + 3df + 

2r2 = 2.4 = cf^ + + a!i^ + + uf + dL'^ 

2r4 = 2.2 = + c?) + + aL'^ + df + M^^^ 

+cf+c^Uaf+af + Mf+d^f^. 



There are also normalization relations among four protons state variables: 



(2.8) 



+ + 4^^W + + 2a^'^ + 2^*^ + 2^^ = 1 

(i = 1 ~ 4), (2.9) 

where the superscript (i) denotes the sublattice number. Finally, by combining eq.(2.5) and eq.(2.7) 
the variational free energy G in the cactus approximation is obtained by 

G = U -TS. (2.10) 

§3. Response to Staggered Electric Field 

Let us calculate the staggered susceptibility to study the properties of phase transition of the 
system. The staggered field is applied so as to induce the anti-ferroelectric order. Since the anti- 
ferroelectric structure is assumed to occur along the a-axis, the electric field E' = —E is applied to 
the sublattice 3 and 4 in addition of E to the sublattice 1 and 2. Since the sublattice 3 and 4 are 
equivalent to 1 and 2 except having a sublattice polarization in the opposite direction, the present 
system is reduced to the one sublattice problem and we have following relations 

,v / _ ^\ / _ ^„ _ ^2 = 1 ~ 4), 

a±, 
d±, 

Pi = P3 = q2 = q4 = n = r2 = S3 = 54 

= co + C2 + ao + a+ + Sd+ + d_, 

P2 = P4 = qi = qs = rs = r4 = Si = S2 

= Co + C2 + ao + a- + d+ + Sd-. (3.1) 





= do, 


% = 


■■ Co, 




— 


(3) 
= «T 


(4) 
= «T 


4" 


= df 
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And the sublattice polarization conjugate to the staggered field is now given by 

p _ p(l) ^ p(2) ^ _p(3) ^ _p(4) 

= a+ - a_ + 2(d+ - d_) = pi - p2. 



The free energy G of the system takes a form: 



A;bT 



^ini- 



-^ + [-2(i±^lnl±^ + ^ 
/cbT ^ ^ 2 2 2 

+(a+ In c-i- + a_ In a_. + 2ao In 

+2co In Co + 2c2 In C2 + 4d+ In (i+ + 4(i_ In d_ )] 

+7^ [1 - (a+ + a_ + 4(d+ + d_) + 2ao + 2co + 2c2)] , 
Kb J 



(3.2) 



(3.3) 



where 7 in the last line is the Lagrange multiplier to make all the state variables 
a+, a_, ao, d_, Co, C2 independent. Under the staggered electric field the thermal equilibrium 
state is obtained from the minimum condition of the free energy: = = = -^j- = = 



dG 



dG 
dc2 



0. The state variables are solved in terms of uq, the electric polarization P and the 



staggered field E as follows: 



Co = Voao, C2 = mao, 
a+ = aoAph^ , a_ = ao{Aph'^)~^ , 
d+ = aor]iAp2 h, d- = aoTjiAp^^ ^ 



(3.4) 



where rjo = exp (— eo/A;BT), r]i = exp (— ei/^bT), rj2 = ex.p {—e2/kBT), h = exp and an 



abbreviation is defined as 



Af 



^ + ^.exp(2DP) {D-^''' 



1-P ' ' k^T' 

Further, oo is determined by a normalization condition as 

2 + 2r?o + 2??2 + Aph^ + Ap'^h'^ 



(3.5) 



ao 



+4r]i{Ap^ h + Ap ^] 



-1 



(3.6) 



Substituting eq.(3.4) and eq.(3.6) into the variational free energy of eq.(3.3), the equilibrium free 
energy Ge under the staggered field is given by 



4an 



A 



Without staggered field this expression has been obtained by Ishibashi.i) The sublattice polarization 
P under the staggered field is obtained as a self-consistent equation from eq.(3.2): 

P = ao \Aph^ - Ap^h-'^ + 2rii{Ap^h - Ap'^h'^) . 

(3.8) 
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The same self-consistent equation for P is also obtained by the thermodynamic relation NfiaP = 
]_aG 

In order to investigate the phase transition properties we expand the above equation (3.8) up to 
the 3-rd order of polarization P and to the linear order of staggered field E in the neighborhood 
of the transition temperature: 

A2{T)P + A,{T)P^ + . . . _ 4(1 + n,)^^ = 0, 

(3.9) 

where A2{T) and A4{T) are given by 

A2{T) = 27]o + 4r/i + 2tj2 - 4(1 + ?7i)D, 



A^iT) = -AD^ -—- + 2riiil + 3D + D^- — 



(3.10) 



From the view point of Landau's phase transition theory the order of the phase transition is classified 
as follows, (i) The phase transition undergoes the second order transition at Tq if ^2 (To) = 
and A4^{Tq) > 0. (ii) The phase transition undergoes the first order transition at Tc{> Tq) if 
^2(^0) = and ^4(To) < 0. The boundary between the first and the second order transition is 
called the tricritical point Tt{= Tq) when A2{Tq) = and A4^{Tq) = 0. The phase diagram in the 
El — T space is shown in Fig. 4. Hereafter the parameters Eq = 6kB,E2 = 1000 A;b, A = 232 A;b are 
taken. It can be seen that the energy parameter Ei representing HPO4 and H3PO4 controls the 
order of phase transition. As Ei/eq is increased, the second order phase transition for small Ei/eq 
changes into the first order transition. In any case the spontaneous sublatticc polarization Pq of 
the anti-ferroelectric state is given by the self-consistent equation (3.8) without electric field: 

Pq = ao[Ap^ - Ap^-^ + 2rii{Ap^^ - Ap^-^)], 

(3.11) 

where clq is a thermal equilibrium value of oq in eq.(3.6) without external field. We solve eq.(3.11) 
numerically and show the spontaneous polarization Pq versus temperature in Fig. 5. On the other 
side, the staggered susceptibility Xs of the system is obtained from eq.(3.8) as the linear response 
AP from Pq induced by the staggered field E 

AP _ Ni^l 1 



E^o E kBT - (^ + D 



l-Pr 







(3.12) 



with 

X = 2do \Apq + Ap^-^ + r]i{ApJ + Ap^~^ 



2Pl 
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Especially, in the paraelectric phase the staggered susceptibility is simplified into 

The staggered susceptibility Xs{T) is shown in Fig. 6 for each case of the second and the first 
order transition. It is noteworthy that when we choose D = 0, the staggered susceptibility never 
diverges. Without a finite dipole-dipole interaction parameter A, we would have the paraelectric 
state down to the zero temperature owing to geometrical frustration structure of the system. 

§4. Response to Uniform Electric Field 

Let us study the response of the present system to a uniform electric field along the a-axis. The 
same external electric field E' = E \s applied to the sublatticc 3 and 4 as well as to the sublattice 1 
and 2. Contrary to the case of staggered electric field, in the anti-ferroelectric state the sublattice 
1 and 2 are non-equivalent to the sublattice 3 and 4 in response to the uniform electric field. The 
present system is reduced to the two sublattice problem and we have following relations: 

4), 



(4.1) 





= do, 


cy = co, cy = c2 (z = 1 ~ 




— a_t ■ 


(3) (4) / 

= a±, = aj_ = a-^, 


4" 


= 4^ ^ 


= d±, = = d':^, 


Pi 


= q2 = 


Co + C2 + ao + a+ + 3d+ + , 


Pi 


= qi = 


Co -1- C2 -1- ao -1- a_ -1- 0?+ -1- 3d-, 


P3 


= 94 = 


Co + C2 + ao + a'.^. + 3d'_,_ -|- d'_, 


P4 


= 93 = 


co + C2 + ao + a'_ + d'_^ + 3d'_, 


2ri 


= 2S3 = 


= 2r2 = 2s4 



= 2co + 2c2 + 2ao + a+ + a'_^_ + 3{d+ + d'^) + (d_ + d'_), 
2rs = 2si = 2r4 = 2s2 

= 2co + 2c2 + 2ao + a- + a'_ + (d+ + d+) + 3{d- + d!_). 

And the sublattice electric polarizations along the a-axis are defined as 

P = p(i) = = _ o_ + 2{d+ - d_) =pi- P2, 
p> ^ _p(3) = _p(4) ^ _ + 2(4 - = P3 - PA. 

(4.2) 

The variational free energy G in the uniform field is rewritten as 

= ^[t-tt; + T-?f{d++d-+d+ + d_) + 



-hl^PP' _ ^"-^ (p _ p'\ 
k^T 2kBT^ ' 
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"^[(Pi liiPi +P2 lnp2) + 2(ri Inri + rs Inra) 
+ (P3 lnp3 +P4 lnp4)] 

+ — [a+ In a+ + a_ In a_ + a'j^ In a'|_ + a'_ In a'_ 
+4ao In oo + 4co In cq + 4c2 In C2 
+4(d+ In d+ + d_ In d_) + 4(d'^ In + d'_ In 
+7^[1 - (a+ + a_ + 4(d+ + d_) + 2ao + 2co + 2c2)] 

+ 7^[1 - (aV + a'_ + 4(dV + d'_) + 2ao + 2co + 2c2)], 



(4.3) 



where 71 and 72 are Lagrange multipliers which are determined by the normalization relations (2.9) 
with the help of eq.(4.1). The thermal equilibrium is determined by the minimum condition of the 
free energy with respect to independent variables: ^ = -§^ = ^ = -§^ = -§^ = -§^ = ^ = 
^F" ~ ^0 ~ ^ ~ ^ ~ ^" '^'^^^^ relations are regarded as a set of equations of state in the present 
system under a homogeneous electric field. In Appendix A it is shown that all the independent 
variables are expressed in terms of sublattice polarizations P and P' under the external electric 
field E. Since the linear response of sublattice polarization to the uniform field E is written as 
AP = P — Pq for the 1 and 2 sublattices and AP' = P' — Pq for the 3 and 4 sublattices, the uniform 
susceptibility Xh is defined as 

Xh = hm —u„ . (4.4) 

Since concrete calculations are complicated and tedious, details of the calculation and the final result 
are given in Appendix A. Here we mention only the homogeneous susceptibility of the system in 
the para-electric phase 

4 + 4r?i 

The numerical results of the uniform susceptibility Xh against temperature are shown in Fig. 7 for 
each case of the second order and the first order transition. 

However, the present uniform susceptibility in the first order transition is completely different 
from experimentally observed one© The observed susceptibility shows the hysteresis phenom- 
ena versus temperature. Until now, we have obtained the susceptibility versus temperature by 
determining the global minimum of the variational free energy at thermal equilibrium for each 
temperature. However, since the heating and cooling process are done at the rate of finite time 
in the experiment, we cannot observe the susceptibility at thermal equilibrium. In order to solve 
this discrepancy, we apply the iteration method called the NIM (natural iteration method)!' for 
numerical calculations of the CVM. The present iteration method leads us to the local minimum of 
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the variational free energy depending upon the initial condition. It should be noted that the local 
minimum is not always the global minimum of the free energy at each temperature. With use of 
the local minimum under the previous temperature as a starting point the sweeping of temperature 
continuously gives the hysteresis curve in the heating and the cooling process (Fig. 8). In order 
to explain the hysteresis phenomena versus temperature we show the temperature change in the 
variational free energy in Fig. 9. Let two temperatures at which a drastic change of uniform suscep- 
tibility occurs in Fig. 8 be defined as Ti and T2 (Ti < < T2). In the heating process, (1) when 
r < Ti, A in Fig. 8 exists in a local minimum denoted by Aq in Fig. 9, (2) when Ti < T < T2, 
B in Fig. 8 exists still in ^0 though another local minimum denoted by Bq appears in Fig. 9, (3) 
at T = Tc the three local minima have the same value and (4) at T = B in Fig. 8 jumps up 
to D in Fig. 8 because the state at the unstable Aq in Fig. 9 tumbles into the stable Bq. In the 
cooling process from D the similar free energy change occurs as a reversible process of the above. 
The actual usage of the NIM is explained more in detail in Appendix B. 

We also calculate the hysteresis curve for the net polarization AP = (P — P')/2 versus homoge- 
neous electric field. The numerical result of AP — E curve is shown in Fig. 10. We can see clearly 
the transition from the anti-ferroelectric state to the polar state at which each polarization in two 
sublattices points to the same direction at E = Ei or E = —Ei. 

§5. A Summary and Discussions 

We applied the cactus approximation of the CVM to the Ishibashi model for the hydro- 
gen bonded ADP-type crystal. The properties of the ADP-type crystal in external electric fields 
were intensively investigated. After re-deriving the variational free energy for the Ishibashi model, 
the equation determining the polarization and the susceptibility in a staggered electric field were 
studied to find the order of the transition. The energy parameter Si characteristic of HPO4 and 
H3PO4 determines the properties of transition, though ADP undergoes the first order paraelectric- 
antiferroelectric phase transition in experiments. We also calculated the susceptibility to a ho- 
mogeneous electric field at thermal equilibrium. The calculated susceptibility does not show any 
hysteresis even in the parameter region of the first order transition, while the homogeneous sus- 
ceptibility observed in experiments shows the hysteresis phenomena in the heating and the cooling 
process. In order to overcome this discrepancy the homogeneous susceptibility in the local minimum 
of free energy was calculated and the result is in qualitatively good agreement with the experiments. 
The hysteresis curve is well explained by utilizing the local minimum in the variational free energy. 
Further, though the typical hysteresis curve of the net polarization versus homogeneous electric 
field have not been found in the ADP experiments to our scarce knowledge, we also calculated the 
hysteresis curve depending upon the external electric field with the same idea. The results are the 
ones expected from the Landau theory of the phase transition in the external field. 
Though in the Ishibashi model the dipole-dipole interaction is included to induce an anti- 
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ferroelectric transition, the present model without dipole-dipole interaction has essentially the 
geometrical frustration and is similar to the spin ice system with ice rule. We will discuss the spin 
ice system from the same view point of the cluster variation method (CVM) in the near future. 

The calculation of the dynamical susceptibility for the ADP-type crystal by the dynamical cluster 
variation method0) is also in progress in order to compare with the experimental data and that0) 
for KDP. 
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Appendix A: Calculation of Uniform Susceptibility 

From the minimum conditions for the free energy G in eq.(4.3) 

dG _ dG _ dG _ dG _ dG _ dG _ dG _ dG _ dG _ OG _ dG _ ct fl^o fr^U^,„inn- 

s^-s^-saT~33r-5^-a^-ad:;:-ad^-^-gs^-^-u, tne toiiowmg 
equations are obtained: 



ao 



[2 + 2t^o + 2m) + tRiR2' 
Co = ??oao, 



C2 = r?2ao 
^1 

«+ 



^^iao(pine2^(^^-P^))/.2^ 
d+ = ^taovi {plp2rfrs) ~' e^'^P'-^'h, (A-1) 



J 



R2 
i?2, 



a'_ 



R2 



R2 



d'_ = ^taovi (p3plrirl)K-''^P^-P'h, 



where t, Ri and i?2 are defined by 

t 



1 ! 

14 



[piP2P3P4.{nr3)'^] 
+ 4r/i {plpirfrs^K^^P^-P^^h-^ 
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(A-2) 



ill 



(A-3) 



A set of self-consistent equation for sublattice polarizations P and P' are obtained from 

P = Pi - P2 = a+ - a- + 2{d+ - d-), 
P' =P3-P4 = o'+ - a'_ + 2(dV - d'_). 

Substitutions of above relations into eq.(A-3) lead us to 

P = ^aot\pirie^''^P^-P^^h^ - p^r3e-^D{p,-p,)^-2 

it2 

+ 2ryi[(p?p2r?r3)^e^(f3-P4);j 
R2 

where it should be noted that Pi,P2,P3,P4:, ''ij ^3 are expressed in the sublattice polarization P and 



(A-4) 



Pi = \i^ + P), P2 = lil-P), 



p, = ^{l + P'), p4 = l(i-P'), 



1,, , P + P\ 



2 
1 

2' 

rs = 77(1 ^ — )• 



(A-5) 



2'* 2 

Thus, eq.(A-4) is the self-consistent equation for P and P'. Since the linear response of sublattice 
polarization to the uniform external field E is defined as AP = P — Po and AP' = P' — Pq, the 
uniform susceptibility Xh is obtained from eq.(A-4) 



as 



where 



^. N ^AP-AP' NfilQi 
^^ = l^oT^« E =k^Q-2 



^ 

dllAp, + Ap^ + 4rjiiAp,-2 +Ap,--2)]' 



(A-6) 
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Appendix B: Natural Iteration Method 

The natural iteration method (NIM)i is one of the method for determining state variables from 
the minimum condition of the variational free energy. Here the NIM is briefly explained from 
starting with eq.(yl-l). It should be noted that bond variables are characteristic of one proton 
configuration while PO4 cluster variables are characteristic of four protons configuration. Four 
protons state variables in eq.(^-l) the left hand side are expressed in terms of only bond state 
variables in the right hand side. Let a temperature T and an external field E with energy parameters 
Eq, El, £2 and A be given. Further, we note that the arbitrary given values of polarization P and 
P' are equivalent to giving the values of bond state variables through eq.(A-5). When the bond 
state variables are substituted into the right hand side of eq.(A-l), the values of four proton state 
variables are naturally calculated. On the contrary, the geometrical relations under the uniform 
electric field 

Pi = 92 = Co + C2 + ao + a+ + 3(i+ + d_ , 
P2 = 91 = Co + C2 + ao + o- + c?+ + 3(i- , 
P3 = 94 = Co + C2 + fflQ + a'j^ + M'j^ + d'_ , 
P4 = 93 = Co + C2 + oq + a'_ + d'_^_ + 3d'_ , 

2ri = 2s3 = 2r2 = 2s4 (B-1) 
= 2co + 2c2 + 2ao + a+ + a'+ + 3(d+ + + (d_ + d'_), 
2r3 = 2si = 2r4 = 2s2 

= 2co + 2c2 + 2ao + a_ + a_ + {d+ + ci+) + 3{d- + d'_) 

give values of bond variables. Thus, these bond variables again determine polarizations P and 
P'. Accordingly, this cycle can be repeated until the convergence within some accuracy is reached. 
Actually it is rigorously proved that as the iteration proceeds, the free energy is always decreased 
toward a local minimum which is not always the global minimum of the free energy. This property 
can be fully utilized in problems of hysteresis phenomena. 
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Fig. 1. The projection of atomie arrangement of ADP-type crystal on (001) plane. The number described in a PO4 
tetrahedron represents the relative height of a PO4 tetrahedron. 
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Fig. 2. Four sublatticies and hydrogen bonds connecting them. Two open circles on a hydrogen bond represent two 
stable points of a proton. 
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Fig. 3. Energy, alloted dipole moment, and occurrence probability of proton configuration around PO4 tetrahedron 
belonging to i-th sublattice. 
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Fig. 4. The phase diagram in ei — T space with eo = 6fcB, £2 = lOOOfee, A = 232fcB. The dotted Une and soUd hne 
represents the first order and the second order phase transition temperature, respectively, and the boundary circle 
represents the tricritical point (TP). 
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Fig. 5. The temperature dependence of antiforroclcctric spontaneous polarization Po- The dotted hne (ei = 
lOfceisecond order) and the sohd hne (ei = 200feB:first order) are shown. 
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Fig. 6. The temperature dependence of staggered susceptibility Xa- The left figure is for the second order case of 
ei = lOfes and the right one for the first order case of ei = 200A;b- 
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Fig. 7. The temperature dependence of uniform susceptibility Xh- The dotted hne (ei = 10fcB:second order) and 
the sohd Hne (ei = 200/cB:first order) are shown. 
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Cooling process: D — ► C — ► A 
Heating process: ^ — ► B — ► D 



Fig. 8. Hysteresis phenomena of uniform susceptibility Xh versus temperature T in the case of si = 200fcB- The soUd 
Une corresponds to coohng process and the dotted Une corresponds to heating process with Ti = 0.9161Tc,T2 = 
l.OOlSTc 
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Fig. 9. Free energy profile. When T < Ti there are two minima at Ao, Ti < T < T2 three minima at Ao, Bo, T2 < T 
one minimum at Bo with ei = 200A;b. 
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